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Abstract The tangent dynamics of the Lyapunov modes and their dy- 
namics as generated numerically - the numerical dynamics - is considered. 
We present a new phcnomcnological description of the numerical dynamical 
structure that accurately reproduces the experimental data for the quasi-one- 
dimensional hard-disk system, and shows that the Lyapunov mode numerical 
dynamics is linear and separate from the rest of the tangent space. Moreover, 
we propose a new, detailed structure for the Lyapunov mode tangent dynam- 
ics, which implies that the Lyapunov modes have well-defined (in)stability 
in either direction of time. We test this tangent dynamics and its derivative 
properties numerically with partial success. The phenomenological descrip- 
tion involves a time-modal linear combination of all other Lyapunov modes 
on the same polarization branch and our proposed Lyapunov mode tangent 
dynamics is based upon the form of the tangent dynamics for the zero modes. 



1 Introduction 

We consider the dynamics of the Lyapunov modes found in two-dimensional 
hard-disk systems. These modes are a subset of the so-called Lyapunov vec- 
tors, which are particular phase space perturbations whose (in)stability is 
characterized by a Lyapunov exponent. The modes themselves exhibit spatial 
and time modality and study of the Lyapunov modes has been consequently 
mainly motivated by their possible connection to hydrodynamic fluctuations 
in fluids. 

The empirical observation of these modes [?,?,?] has led to various stud- 
ies of their physical origins based on several different approaches. Firstly, 
random matrix approximations have demonstrated the importance of trans- 
lational invariance in the phase space [?], while a combined random matrix 
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and master equation approach [?] has predicted the existence of the Lyapunov 
steps, which are the sets of degenerate Lyapunov exponents characterizing 
the Lyapunov mode (in)stability, and of the spatial modality. Alternatively, 
application of kinetic theory in Ref. [?] exhibited a connection between mode 
generation and conservation laws, and predicted polarization of the modes 
into longitudinal and transverse polarization branches. Further, describing 
Lyapunov modes as Goldstone modes [?] arising from symmetry breaking 
produced approximations to the wavenumber dependence of the Lyapunov 
exponents and the form of the modes. 

The abovementioned approaches focus mainly on the generation and spa- 
tial structure of the Lyapunov modes. The underlying general theory of the 
tangent dynamics of the Lyapunov modes is well-known (see e.g. [?] or [?] 
for a review), while study of the dynamics of the Lyapunov modes as they 
are generated numerically - the numerical dynamics - has included: observa- 
tion of time-periodicity in the longitudinal modes [?,?,?,?]; prediction of the 
consequent mode propagation speeds for low densities (see e.g. Refs. [?,?]); 
and has revealed a connection between the time oscillating period and the 
momentum autocorrelation function [?]. Refs. [?,?] have considered the nu- 
merical dynamics of longitudinal modes in terms of rotations in the tangent 
space, but the numerical dynamical structure of the modes has not been ex- 
amined in the same level of detail as the spatial structure or the degeneracy 
of the Lyapunov exponents. 

In the present work, we present a phenomenological description of the 
numerical dynamical structure of modes in the quasi-one-dimensional system, 
which accurately reproduces the numerical dynamics of its Lyapunov modes. 
This description demonstrates that the mode numerical dynamics is linear 
and appears to separate from the numerical dynamics of the other, non-modal 
Lyapunov vectors. Moreover, we propose for the first time a detailed structure 
of the mode tangent dynamics. This mode tangent dynamics possesses a 
symplecticity property, which is shown to imply that the Lyapunov modes 
lie in the so-called Oseledec decomposition spaces, and so have well-defined 
(in)stability in either direction of time. We test this tangent dynamics and its 
derivative properties numerically with partial success. The essential feature 
of the phenomenological description is that the numerical dynamics of a 
Lyapunov mode is a time-modal linear combination of all other modes on 
the same polarization branch. Our mode tangent dynamics is based on this 
phenomenological description and an extension of the zero mode tangent 
dynamics. 

2 Lyapunov Vectors and Modes 

2.1 Tangent Space Dynamics 

Consider an r-dimensional system of TV identical particles, with phase space 
M ~ M 2rAr . For a state £ <E M the dynamics £(£) is represented by the phase 
flow </>* : M -> M. That is, = where 



i = (p,q) = (Pi, ■ ■ ■ ,PN,qi, ■ ■ ■ ,qw) 



(2.1) 
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and pi (qi) is the momentum (position) of the ith particle. 

In the present work we consider a system of hard-spheres, which interact 
only by instantaneous elastic collisions. There is no external potential in this 
system and we neglect tangent collisions so that this phase flow is differen- 
tiable. A small perturbation St; of the state £ then has dynamics (to first 
order) 

5£(t) = Dtf&t, (2.2) 

where D^ft = deft/dt; is the 2rN x 2rN matrix of partial derivatives of the 
phase flow. The map D^cj) 1 is called the tangent flow, and we may consider 
£>£</>' : T ? M -> IWjM, where T$M is the tangent space to M at £. The 
dynamics in Eq. (|2.2[) therefore implies a perturbation 6£ may be considered 
to be a tangent vector, which is an element of the tangent space. We hence- 
forth adopt the convention of calling 5£ a tangent vector. The components 
of 5£ are written 

5£ = (Sp,Sq) = [Spi, . . .,Sp N ,Sqi, . . .,Sq N ) . (2.3) 



2.2 Lyapunov Exponents and Vectors 

A detailed summary of the following theory of Lyapunov vectors and the 
numerical scheme is presented in Appendix [Al Firstly, it is well-known that 
one may choose an orthogonal basis of the tangent space such that in the 
positive time limit t — > oo, each basis element 6£t satisfies 

\\5Z+(t)\\^ex P (X^k)\\5$\\ , (2.4) 

for k = 1, ... , 2rN and 1 < j k < I (see App. lA.lj) . Here I < 2rN is the 
number of distinct exponents A which are called Lyapunov exponents 
and are ordered such that A^ > . . . > \"h The ordered set of Lyapunov 
exponents is called the Lyapunov spectrum, and we say the tangent vector 
is 'associated with' the Lyapunov exponent X^ k \ 
The notation \W here denotes a distinct Lyapunov exponent. A par- 
ticular exponent A^ may have multiplicity m(j), such that jh = ■■■ = 
ik+m{j) = j- The associated basis elements ■ ■ ■ > ^^t+mtj) then span a 
subspace of TjM, whose elements are themselves all associated with A^. 
There are therefore many possible (pairwise linearly independent) choices of 

5^,..., 5£+ +m[j) for m(j) > 2. We label this subspace C/^j C T$M, and call 

the elements of + Lyapunov vectors. 

The system under consideration is Hamiltonian, so that the dynamics (j) 1 
is reversible and we may similarly consider the behavior of tangent vectors in 
the negative time limit, t — > —oo (see Sec. IA.2]) . In this case, the Lyapunov 
spectrum is -A« > ...-AW, and we may choose another orthogonal basis 
of T e M such that 

||^WII-exp(-A^)||^7|| . (2.5) 
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There similarly exist subspaces _ C T^M whose elements are associated 

to — X^i for t — > — oo. 

Notably, a Lyapunov vector associated to — \^ for t — ► — oo is not nec- 
essarily associated to A^ for t — ► +00, so that ^ f/i^i m general. 
However, it is possible to decompose the tangent space into subspaces whose 
elements do have well-defined Lyapunov exponents in either direction of time 
(see App. IA.3p . These subspaces are called Oseledec spaces, denoted 
and a tangent vector 5£ £ satisfies both Eqs. (|2.4| and (|2.5[) with ex- 
ponent X( j \ The Oseledec spaces are also covariant: they are preserved by 
the tangent dynamics, so that Dg0*Wt = wit'Ly 



2.3 Numerical Scheme 

A set of 2rN orthonormal Lyapunov vectors may be generated numerically by 
the numerical scheme of Benettin [?] and Shimada [?] , which is used primar- 
ily to generate the Lyapunov spectrum. This scheme (see App. IA.4[) involves 
evolving a set of randomly chosen, linearly independent tangent vectors for- 
wards in time under the tangent flow, whilst periodically performing Gram- 
Schmidt orthonormalization. The resulting set converges to an orthonormal 
Lyapunov vector basis. Importantly, the Lyapunov vectors generated by this 
scheme are strictly associated with —A"' for t — > — oo, rather than with +A^' 
for t — > +00 as would be expected naively. 

The Hamiltonian nature of the system further guarantees (see App. IA.5|) 
that there is always a Lyapunov exponent equal to zero (that has even mul- 
tiplicity) and all other distinct exponents are paired signwise, so 

A a) = _ A a-j+D. (2.6) 

It follows that the number of distinct exponents I is always odd and that 

A ([/+l]/2) _ Q 

Furthermore, a basis of 2rN Lyapunov vectors may be uniquely specified 
by considering only the first rN Lyapunov vectors generated by the numerical 
scheme (see App. lA.5"|) . These rN Lyapunov vectors span U±^(B- • ■®U^ l ^ i ^ 2 ' 

together with half of U^ 1 ^ 1 ^ 2 ^. We henceforth call these Lyapunov vectors 
the positive Lyapunov vectors, since they are associated with exponents A > 
0. The remaining rN Lyapunov vectors associated with A < are called 
negative Lyapunov vectors. 



2.4 Quasi-One-Dimensional System 

Henceforth we consider only two-dimensional rectangular systems of N parti- 
cles with rectangular boundaries, so in Eq. (|2.1[) the components Pi,qi £ K 2 . 
We use Cartesian coordinates (x, y) in this system and let q = (x,y) and 
P = (Px,Py)- For simplicity, we will mainly consider a special case of the rect- 
angular system called the quasi-one-dimensional system. This system is a 
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two-dimensional system which is sufficiently narrow so that particles cannot 
pass one another (see Fig. [T]). 




X 



Fig. 1 Schematic of an N hard-disk quasi-one-dimensional system. The height 
L y is sufficiently small that the disks cannot pass one another. We choose the 
coordinate origin to be located at the bottom left corner of the system, and system 
boundaries are denoted by dashed lines. The pairs of opposite boundaries x = 0, L x 
or y — 0, L y may possess either periodic (P) or hard (H) boundary conditions. We 
call the system with hard x = 0, L x boundaries and periodic y = 0, L y boundaries 
the HP quasi-one-dimensional system. 



2.5 Zero Modes 

In the case of periodic boundary conditions, the rectangular system has full 
translational spatial and time symmetry. There are then (see App. IA.6[) six 
Lyapunov vectors associated with Ad' +1 l/ 2 ) = in either time limit. The 
first of these Lyapunov vectors corresponds to uniform translation in the x 
direction and is <J£i = l/A^ 1 / 2 (0, 0, ei, . . . ,ei), where e\ = (1,0) is repeated 
N times and is the N dimensional zero vector. We simplify this notation 
by writing (see Eq. 12.3ft 

5£= (•••,(<5p»*.<*P» tf ),- ••,•••,(%„)%„),•••) = {Sp x ,Sp y ,Sq x ,Sq y ) , 

(2.7) 

so 8£i — l/A^ 1 / 2 (0,0, 1,0). The other five Lyapunov vectors are then <5^2 = 
1/^/2(0,0,0,1), *£a = l/|b||(0,0,P*,Ps,), = 1/^V 1/2 (1, 0,0,0), % - 
l/7V 1 / 2 (0, 1, 0, 0) and S^e = l/l|p||(Px,P y , 0, 0), which correspond respectively 
to translation in the y direction, time translation, a uniform momentum boost 
in the x direction, a uniform momentum boost in the y direction, and an 
energy shift. Note that in the case of e.g. $£3, putting Sp x — p x means that 
dpi* = Pi x ■ These Lyapunov vectors are collectively called the zero modes. 
One finds that linear combinations of 8£x, 6^,2 and 6^3 (££4, <5£s and S^e) 
belong to the positive (negative) Lyapunov vectors (see Sec. 12. 3[) . 

In the case that boundary conditions at one pair of opposite boundaries of 
the rectangular system are reflecting rather than periodic, then one symmetry 
is broken and the number of zero modes is reduced by two. For a quasi-one- 
dimensional system with mixed period and reflecting boundary conditions, 
we thus expect four zero modes. We assume in the present work that time 
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symmetry is never broken and total energy is always conserved, so there are 
always at least two zero modes for the system. 

2.6 Lyapunov Modes in the Quasi-One-Dimensional System 

Consider a positive Lyapunov vector, <5£ € T^M of unit norm at some fixed 
time t. Generally in any r dimensional system, the Lyapunov vectors gen- 
erated by the numerical scheme are strongly dynamically localized (see e.g. 
Ref. [?]). This means that for a small, rapidly changing index subset of the 
integers 1, . . . , 2rN the components Spi,5qi ~ 1 and 5pi, 5qi <C 1 elsewhere. 

However, Lyapunov exponents close to zero exhibit a degenerate step- 
like structure, and the Lyapunov vectors associated with these exponents 
are found to be delocalized over the N particles [?,?,?]. More precisely, the 
components 5pi (Sqi) are approximately linear combinations of 4> p (pi) (<fi q {qi)) 
and Piipp(qi) (Pi^> 9 (<?»)), where <f> p ,<fi q ■ K 2 — > M 2 and ip p ,tp q : M 2 — > K are 
modal functions. The functions 4> p or 4> q are called modal vector fields and 
Piipp or Piipq are called a scalar modulation of the momentum field. The 
Lyapunov vectors with this modal property are called Lyapunov modes, and 
in this work will be denoted by Srj. The Lyapunov modes may be represented 
by the fields <p p , <f> q , tp p and tp q which approximate them. We call this the 
vector field representation of the Lyapunov mode. 

We call the modes belonging to the positive Lyapunov vectors the positive 
modes. For a general two dimensional system, the zero modes belonging to 
the positive modes - the positive zero modes - are linear combinations of 

»*-(J), **-(»). *» ^ 

Note that we write the positive zero modes in the vector field representation, 
and we only write the spatial components, as the momentum components 
are all zero. 

For a general two dimensional system, the modal scalar or vector fields 
that well-approximate the non-zero Lyapunov modes have a certain number 
of nodes in the x or y direction. We define n x (n y ) to be the number of nodes 
in the x (y) direction of this Lyapunov mode and define the wavevector of 
the mode to be 

k= (k x ,k y ) = (nn x /L x ,Tm y /Ly^j , (2.9) 

where L x and L y are the spatial sizes of the system. There may clearly be 
several modes with the same wavenumber ||fc||, but for a particular wavenum- 
bcr, one finds that the corresponding positive modes are associated with 
Lyapunov exponent 

AwCrllfeH, or X^C L \\k\\ , (2.10) 

where Ct and Cl are some scalar constants of the system (see e.g [?] or [?] ) . 
The modes associated with Ci||fc|| are called longitudinal modes, while those 
associated with Ct 1 1 1 1 are called transverse modes. Generally there are twice 
as many longitudinal modes as transverse modes for any particular ||fe||: a 
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longitudinal mode may be written as a5r]L + PSrjp, where 8t]l is formed 
by a modal vector field called an L-mode and 5t]p is a formed by a scalar 
modulation of the momentum field called a P-modc; while a transverse mode 
is written as St/t, a modal vector field called a T-mode. For this reason, 
the longitudinal (transverse) modes are often referred to as LP-modes (T- 
modes) , and the two linearly independent orthogonal linear combinations of 
a particular Stjl and Srjp are called an LP mode pair. 

One notable empirical feature of the non-zero positive or negative Lya- 
punov modes is that the momentum component dp of a mode is found to be 
proportional to the spatial component 8q. That is, we may write 

S V ^(CSq,Sq) , (2.11) 

where C is a scalar. This relation means that we need only specify the spatial 
part of the Lyapunov mode in order to describe the entire mode (see e.g. [?]), 
and henceforth we consider only the spatial parts of the mode unless stated 
otherwise. Equation (|2.1ip also permits us to visually represent a Lyapunov 
mode by drawing each (2-dimensional) Sqi at the corresponding coordinates 
qi. For a quasi-one-dimensional system we may plot instead (8xj) or (Syj) 
against (xj) — jL x /N (see Ref. [?]). 

For convenience, we henceforth set L x = 2tt. For the HP quasi-one- 
dimensional system, which is of particular interest in the present work, the 
vector field representations of the positive T, L and P modes respectively are 
found to be approximately of the form [?] 



SrjT{k x ) = V% kJ cos(k x x) , 
SriL(k x ) = V% \t\ sm(k x x) , 
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Svp(k x ) = if,, I P * ) cos(k x x) , (2.12) 



Ibll \Pv 

for k x = 1/2,1,3/2,.... 



2.7 Inner Product and Basis 

The Lyapunov modes generated by the numerical scheme are orthonormal. 
In the vector field representation, the inner product, which defines orthonor- 
mality, is simply 

(5rii,8rij) = — — / / • 8r]jdxdy , (2.13) 

L X -Ly JQ JO 

where Sr/i ■ 8rjj means a dot product of Srji and Srjj in the vector field repre- 
sentation. Note that by conservation of momentum, under this inner product 

I I p x dxdy = I / pydxdy = , ( 2 -14) 
Jo Jo Jo Jo 
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and the T, L, P and zero modes in Eqs. (12.81) and (|2.12jl have unit norm 
under this inner product. 

/From Eq. (|2. 12[) it follows clearly that (£?7t, <%l) = 0. Further, as ex- 
pected (Srj a (i) , Sr] a {j)) — 0, for a — T, L or P and half-integers or integers 
i j. However, 8rjp(k x ) is not orthogonal to Sr)T(k x ) in general, since 

{Sr] T (k x ),Sr]p(k x )) = — — / p cos(2k x x)dx . (2.15) 
47r l|Pll Jo 

Similarly the zero modes 5rj\ or $772 are not strictly orthogonal to the P- 
modes. Of course, the T (zero) and P (P) modes must always be orthogonal, 
as they are Lyapunov vectors associated to different Lyapunov exponents. 
This difficulty may be resolved by presuming that the p y scalar field is random 
white noise, so that the integral is zero. 

Note also the L and P modes Sr]L(k x ) and 5rip(k x ) are not strictly or- 
thogonal, since 

(5ilL(k x ),5r]p(k x )) = / p x sm(2k x x)dx . (2.16) 

47r lbll Jo 

This may be easily resolved by performing an orthonormalization procedure 
on the L and P modes, which are associated with the same Lyapunov ex- 
ponent and need not be strictly orthogonal. However, in line with our pre- 
sumption for the p y , we instead presume the p x are also white noise, so the 
integral is zero. 

Using notation in Eqs. (|2.8p and (|2.12p . a basis of positive Lyapunov 
modes for the HP quasi-one-dimensional system (see Fig. [1]) is then the set 

Sm^m,Svr(^,Sr]L(^,Sr]p(^j,. . . ' <S? ?T(y),^L(y),^p(y) 

where K = 3m +2 is the number of Lyapunov modes in the system. Note that 
the Srji zero mode is absent, since translational symmetry in the x direction 
is broken by the hard boundary conditions. 



3 Lyapunov Mode Phenomenology 

3.1 Numerical Dynamics 

Let fl l be the matrix whose columns are the AN orthonormal Lyapunov 
vectors <5£* generated at a particular time t by the numerical scheme, with 
J? = Q. We call the dynamics of a Lyapunov vector or mode under the 
numerical scheme the numerical dynamics. This dynamics is not equivalent to 
the actual tangent space dynamics (12. 2|) , since the numerical scheme involves 
both evolution under the tangent flow combined with regular application of 
the Gram-Schmidt orthnormalization procedure. 

The numerical dynamics of the Lyapunov vectors may then be written as 



(3.1) 
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where Nq : T^M — ► T^t^M is the numerical dynamics operator, which 
preserves orthonormality. The notation in Eq. (13. ip indicates operation of N n 
column-wise on the matrix 12, such that if Q = (. . . \Sujj {■•■), then N^Sujj) — 
Suij. Although Nq(Q) is a matrix under this notation, the operator Nq is 
not generally a linear map, as its form is dependent on the choice of J? 
and in particular the ordering of its columns, due to the iterative nature of 
the Gram-Schmidt procedure. We indicate this sensitivity to the choice of 
Lyapunov vector basis and its ordering by the subscript fl. Note that both 
TjM and ?> (?) M ~ R 4N , so that we may consider Nq : R 4N -> R 4N , or 
Nq : R 2 — > R 2 in the vector field representation. 

The Lyapunov modes generated by the numerical scheme may exhibit 
time-periodic behavior. That is, for a Lyapunov mode 8rf at some time t, 
there exists some period r such that 8rf = Sr] T+t . The longitudinal modes 
generally exhibit the most obvious time-periodicity, as can be seen for the 
HP quasi-one-dimensional system in Figs. 4 and 5 of [?]. 

Many previous works (e.g. [?,?]) have noted that the numerical dynam- 
ics of an LP mode pair are well-approximated by a time-modulated linear 
combinations of their corresponding L and P modes at some fixed time, say 
t,Q = 0. For example, putting r = 2n, we might have 

Srf = f^Q » V2cos(i) Q sm(x/2) + (i 

as an approximation to longitudinal mode 197 shown in Figs. 4 and 5 of 
Rcf. [?]. 

Figure [2] shows the contour plots for the predicted LP mode numerical 
dynamics according to Eq. p.2|) . Note that we neglect terms involving p x in 
the x-component of the mode, since we assumed this momentum component 
was random white noise. Similarly, the y-component of an LP mode includes 
only terms with p y coefficients, so we plot instead 8rf y /p y . We will plot all 
other modes in this way. 

The empirical contour plots in Fig. 5 of Ref. [?] reveal notable features 
of the LP mode dynamics such as triangular-like half-height contours that 
have a dimple on their short side. These are not reproduced in Fig. [2] here. 
Hence, at least in the quasi-one-dimensional case, it is clear that the time- 
periodic behavior of the LP mode is not fully described by Eq. (|3.2|) . There 
are therefore two questions which arise concerning the numerical dynamics: 
can a better description of the numerical dynamics be found; and moreover 
is the time periodicity merely an artifact of the numerical scheme or a more 
fundamental property of the tangent dynamics. 

3.2 Separation of the Numerical Dynamics 

The numerical dynamics operator Nq is bijective, so that at some time t 
a Lyapunov vector may be described as some linear combination of all the 
other Lyapunov vectors at some initial time t — 0. However, the partially 
successful description (|3.2j) of an LP mode at any time t involves a linear 
combination of only its corresponding L and P modes at t — 0. Validity of 



J cos(a;/2) , (3.2) 
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Fig. 2 Contour plots of the predicted LP mode numerical dynamics in the HP 
quasi-one-dimensional system according to Eq. (|3.2[) for (a) Sr/l and (b) 5rf y /p y . 
Three contours are shown in each plot: a contour at the half maximum height level, 
1/V2 (dashed line); a contour at the zero level (solid line); and a contour at the 
half minimum height level, — l/v2 (dotted line). 



Eq. (|3.2p implies then that the numerical dynamics of this LP mode pair 
separates from the dynamics of other Lyapunov modes or vectors, and that 
the numerical dynamics is a linear map. We now use this general principle 
to obtain a better description of the numerical dynamics with the following 
claim. 

Let {Srjk}, k = 1, . . . , K be a basis of positive Lyapunov modes ordered as 
in Eq. (|2. 1 7|) . Let £"* be the 2 x K matrix with orthonormal columns whose 
fcth column is the positive mode bify. in the vector field representation. For 
the HP quasi-one-dimensional system, by Eq. (|2.17|) we set the vector field 
representation form for S° = E to be the 2x11 matrix 



5t] T 



8r]L 



^ P (l)\--]s VT Q\8 VL Q\Sr, P Q 



(3-3) 



Claim After some time interval t, the positive Lyapunov modes 8rf k are time- 
modulated linear combinations of only the positive initial modes Sr]k such 
that 

S 1 = SA(t) , (3.4) 

where Ait) is a modal, orthogonal K xK matrix. We call A(t) the admixture 
matrix since it represents linear combinations of Lyapunov modes. We further 
claim the longitudinal (transverse) modes are linear combinations only of the 
initial longitudinal (transverse) modes. 

According to Eq. 13. 4[ the j th column of Ait) , a,j it) , is matrix- multiplied 
with the matrix S to produce the jth mode at any time t. We say that a jit) 
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generates the jth mode, and note that A(t) must be orthogonal since the 
numerical scheme produces orthonormal sets of Lyapunov modes: both S 
and S have orthonormal columns, and hence A(t) must be orthogonal. The 
claim (|3.4[) implies that the numerical dynamics of either the longitudinal or 
transverse modes separates from that of the rest of the tangent space. 



3.3 Admixture Matrix 

In order to test the claim (|3.4p , we now seek to construct the matrix Ait) 
for the HP quasi-one-dimensional system (see Fig. [TJ. For simplicity, we fix 
the period r = 27r, and we only construct explicitly the columns of Ait), 
which generate the zero modes and those modes corresponding to smallest 
wavenumber k x — 1/2. That is we consider only the first five columns of 
Ait). 

The positive zero modes are static under the numerical dynamics (see 
App. IA.6j) . and assuming the transverse modes are also static^], from Eq. (|3.2j) 
an initial guess at the matrix Ait) is the block diagonal matrix 



Ait)= (ai(i)|...|a u (*)) = 



/10 
1 









cos(i) sin(t) 
— sin(t) cos(i) 
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(3.5) 



where the submatrices on the diagonal are rotation matrices with period 2n. 

Now, the desired finer details of the LP modes, such as dimples in the 
half-height contours, appear to occur at double the frequency of the actual 
mode, so we presume then that there are terms such as cos(t) cos(2i) in 
the admixture matrix: we expect a modulation on top of the modulation 
described in Eq. (|3.2|) , so that the former has double the latter's frequency. 
Since A(t) is orthogonal, its columns must have unit norm, so if cos(i) cos(2t) 
occurs in a column then so must cos(i) sin(2i), sin(t) cos(2i) and sin(i) sin(2t). 
Treating Eq. (13. 5p as a starting approximation, and enforcing orthonormality, 
a possible choice of the fourth and fifth columns of Ait) is 



a 4 (t) 



a 5 (t) 



7T 



0. 



7T 



L=(o,o, 

+ a 2 V 

c(i), -s(t), 0, ac(t)c(2t), -as(t)s(2t), 0, as(f)c(2t), -ac{t)s{2t)j 

(3.6) 

= (o,o,o, 

4- rv 2 V 



sit), c(t),0, as(t)c(2t), -ac(*)s(2t), 0, -ac(t)c(2t),as(i)s{2t) 



(3.7) 



1 Fourier analysis of the empirical transverse mode numerical dynamics is re- 
quired to confirm this. We do not address this issue in the present work. 
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where c(t) = cos(t), s(t) = sin(t), and a > is some small real number. 

Plots of the predicted modes generated by these columns are shown in Fig. 
[3l As can be seen by comparison with Fig. 5 in Ref. [?] , the finer structure of 
the mode dynamics, such as the triangular positive and negative half-height 
contours with dimples, appear to be properly reproduced and oriented by 
Eqs. (|3.6|) and (13. 7|) . Further, the zero level contour of the Sr] y /p y components 
appears to match closely the numerical data, although this is more difficult to 
determine due to numerical noise. Equations (|3.6p and (|3.7p therefore appear 
to better describe the numerical dynamics than Eq. (|3.2j) . 



4 Lyapunov Mode Dynamics 

4.1 Full Numerical Dynamics of the Lyapunov Modes 

For a basis of positive Lyapunov modes {Sr/k}, k = 1, . . . , K as in Eq. (|2.17p . 
we may define a corresponding basis of negative modes (see App. IA.5|) by 
putting 

Srik+K = JmSr/k ■ (4.1) 
Here for some integer L, J 2 l is the 2L x 2L matrix 

J2L = (°j "/) , (4.2) 

in which (I) is the L x L null (identity) matrix. 

We define the matrix S = (. . . \Sr]k\ ■ ■ ■) similarly to Eq. (|3.3jl . and then 
by Eq. (|4.ip . (E\J iN EE) = (£?7fe)fc=i ^ s a ma t r i x whose columns form an 
orthonormal basis of all the Lyapunov modes. Note that we consider here 
Sr]k £ K 4Ar , rather than the vector field representation, and the inner product 
of modes is then simply the dot product Srji ■ 5rjj = Sr/fSr/j, rather than an 
integral (|2.13p . The following arguments will, of course, hold for the vector 
field representation and its inner product. Note also that the orthonormality 
of the columns of (S\ J^S) is due to the property 

~ T J iN ~ = 0. (4.3) 

This depends on the fact that positive mode Srjk £ U^_, for some j < 

(I + l)/2. However, JSrjk £ JU^l = J7^I J+1) (see App. EH) is a negative 

mode and the spaces are all orthogonal, so Sr)i ■ JSrjk = for any i and 
k. 

The claim (|3.4p describes the numerical dynamics of only the positive 



modes. However, the numerical dynamics Nq (13. ip is a symplectic operator 
(see Appendix IB"]) such that (J^nS) 1 — J 4 at(S')'. Hence by the claim (|3.4p . 
the negative modes must have dynamics (J^lE;)' = J4jv(H) = J^N^A{i). 
The claim (|3.4p therefore implies that the full numerical dynamics of all the 
modes is simply 

(h|j 4J vS)* = (s\j 4N s)Af(t), where JV(t) = ( A ® J>\ {AA) 
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Fig. 3 Contour plots of the predicted first and second LP mode num erical dyn am- 
ics in the HP quasi-one-dimensional system according to Eqs. (|3.6[1 and (|3.7|) for 
(a) 5rj x of the first LP mode, (b) Sr/ X of the second LP mode, (c) Sr/y/py of the 
first LP mode, and (d) Sri y /p y of the second LP mode. The parameter a = 0.2. 
Three contours are shown in each plot: a contour at the half maximum height level, 
l/y/2 (dashed line); a contour at the zero level (solid line); and a contour at the 
half minimum height level, —l/y/2 (dotted line). Compare this figure with Fig. 5 
of Ref. [?]. 



Since, according to the claim (|3.4[) . the matrix A{t) is orthogonal, then 
so is the 2K x 2K matrix Af(t) . The orthogonality of A(t) also means that 
Af(t) is a symplectic matrix. That is, 

N{t) T J 2K N{t) = J 2K , (4.5) 

as expected (see App. |B|. 
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Further, by Eq. (13. 1| . the numerical dynamics of the modes may be gen- 
erally written as 











[(* 


Jan — ^ 



(4.6) 



where here iV^ 
not matter as 



p4JV 



p4W 



Note that the ordering of [S\J^S) does 
is a column- wise operator (see Eq. I3.1|) . Let TgJf = 

implies that T^AT is an 



spanj^ryfc}^^. Then the numerical dynamics (|4.4 
invariant subspace of Nq, that is N^lT^X) — T^X, or in other words, the 
numerical dynamics on T^X separates from the dynamics on the rest of the 
tangent space. Furthermore, comparing Eq. (|4.4p with Eq. (14.61) it follows 
that iV&[(S| Jw£)] = (3\J4jfE)Af(t), and hence 



Jan- ) 5% 



(4.7) 



Since Sr/k is a basis of T^A", then it must be that the numerical dynamics Nq 
is a linear operator on T^X, provided the claim (|3.4|) is true. 



4.2 Connection to the Tangent Dynamics 



Consider the usual tangent dynamics (|2.2[) after normalization, which we 
henceforth call the normalized tangent dynamics. In general, the normalized 
tangent dynamics of the modes is a linear map and may be written as 



D(t)(S 



,h N S ) = D. 



Jan- diag 



1 



WD^SmW J k=i 



2K 



(4.8) 



We define Srj k (t) = D(t)6r) k . 

Suppose the normalized tangent dynamics on T^X separates from the 
other dynamics, so that D{t)8r] k = [E\JAN3)T>{t){E\JAN3) T 8r]k for some 
matrix T>(t) similarly to Eq. (|4.7p . Then, since the Srjk form a basis of T^X, 
it follows that we have 



noting orthonormality of the columns of [B, \ Jm^z) means 



Jan— ) — J. 



t2K 



(4.9) 



(4.10) 



where I^k is the 2K x 2K identity matrix 

The normalized tangent dynamics of the zero modes clearly separates 
from that of the rest of the tangent space (see App. IA.6]) . From Eqs. (|A.28|) . 
(|A.29|> and (|A.30|) we have that 



D 



2s+2 



Sms+2)v {t) 



(4.11) 
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where Sr/i , . . . , 8rj 2a+2 are the zero modes for a system with s spatial symme- 
_ [I s+1 (3tl s+1 [{(3tf + 1]-V2\ 

We may trivially decompose T>o(t) as follows: 

PoW-/ 2s+ 2^ Q J s+i p t) 2 + 1] -l/2 J . (4-13) 

noting that the identity I2.S+2 is simply the numerical dynamics of the zero 
modes Ao(i). Since the Lyapunov modes appear to be modulations of the 
zero modes [?], we then extend this dynamical form and presume that the 
normalized tangent dynamics of the modes is 

m _ (A{t) \ (Ik ptI K mf + 
V ®-[ A(t))[o lKmr + 1] -i/2 ) , (4-14) 

where we have replaced the numerical dynamics of the zero modes Afo (t) by 
our claimed numerical dynamics AT(t) (|4.4|) for the modes, and adjusted the 
dimensionality of the matrix on the far right accordingly. That is we claim 

_ (A{t) ptA{t)[{ptf + 

is a general form for the normalized tangent dynamics, where A(t) is some 
periodic, orthogonal matrix. Note that Eqs. (|4. 1 3[1 and (|4.14[) are the QR 
factorizations of T>o(t) and T>(t) respectively. This is expected since the nu- 
merical scheme differs from the normalized tangent dynamics only by the 
application of the Gram-Schmidt orthonormalization procedure. 



4.3 Symplecticity of the Normalized Tangent Dynamics 



We may further re- write Eq. ([4TT4]) as V(t) = Af(t)K(t), defining TZ(t) in the 
obvious way. One finds that the matrix lZ(t) is 'almost symplectic', in the 
sense that TZ(t) T J 2K TZ{t) = J 2K [(/3t) 2 + l]" 1 / 2 . By Eq. gl)]) it follows that 
T>(t) is similarly almost symplectic, i.e. 



V{t) T J 2K V{t) = J 2 K[((3t) 2 + 1P 1/2 
Now, the matrix (S\J4^S) has the property that 

T 



J: 



2K, 



and from Eqs. (|4~8|) . (|4~9|) and (j4. 10|) we have explictly that 



V(t) = (s\j m Sy D^fs 



Jan~ )dia; 



\D t ^Sxik\ 



2 A' 



fc=l 



(4.16) 



(4.17) 



(4.18) 
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Separability of the normalized tangent dynamics (|4.9jl implies that its projec- 
tion onto X, which is (S\ J4nE)(E\ J^m^) T D^(f) t {S\ J^S), is simply equal 
to D i <j) t {S\J iN S). Hence from this together with Eqs. (|4. 18[) and (|4.17p we 
have that Eq. (|4.16[) is equivalent to 

f 1 ~\ 2K f 1 V K 

(4.19) 

which in turn holds if and only if 

logHA^feH +log||£) e ^J4 J v%|| = log{[(/3t) 2 + 1] 1/2 } (4.20) 
for k= 1,...,K. 



4.4 Lyapunov Mode Stability 

The empirical multiplicity of the zero modes implies (see App. IA.7[) that the 
spaces of Lyapunov zero modes satisfy 

where m = (I + l)/2 and is the zero Oseledec space. 

Suppose that the symplecticity property (|4.20[) of the Lyapunov modes 
holds. We call the spaces U^_, whose elements are Lyapunov modes, mode 

spaces, and the symplecticity property then implies that for all <5£ € uj?}_, a 
mode space, 

Urn 1||£>^||+ lim h\D^J5^\ = (4.22) 

since log[(/3i) 2 + l]/|t| -t as ^ ±oo. Certainly Eq. (|4~22|) is true for 
the limit t — > — oo by the general theory of tangent vectors, as J?7fc is the 
Lyapunov vector associated to Lyapunov exponent A& and dr]k+K to — A& 
in the negative time limit (see Apps. IA.ll IA.2[) . However, this result is not 
necessarily true for t — ► +oo, and thus is a particular property of the form 
of our claimed normalized tangent dynamics (|4.15|) . We now show that Eq. 
1)4.220 means that 

U^ = Wj?\ (4.23) 

where j is such that is a mode space. 

In order to show uj?}_ — , it is sufficient to show C Wy'' because 

dimW $ (m) = dimL^™ (see App.[OJ). First of all, it follows from Eq. (|A~T6)) 
that 

Ug>_ C e T e M : t hrn^ i log ||£> S #^|| > } U {0} . (4.24) 
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We then have that, by Eqs. (|A.25|) and (|A.26j) . 



u® = julz 3+1) 



C G T e M : lim 1 log JS£\\ > -X®) U {0} 

= keT 5 M: lim llogllD^'^H < A«lu{0} , (4.25) 

the last step of which follows from the symplecticity relation (|4.22[) . Com- 
paring Eqs. (|4.24|) and (|4.25[) produces 



Ugl C |^ G T 5 M : t Urn^ 1 log ||£>^e|| = A^ J U {0} 



(4.26) 



for a mode space U^_. But by Eqs. (fA~TT|l and jll2| , for any eigenspace it 
is the case that 



C j# G T^M : k^i log 11^*^11 = -A (J) j U {0} , (4.27) 



and thus 

t^l C € T^M : t lirn^ ^ log ||£> € 0'^|| = ±A^ Ju{0} = . (4.28) 

Note that the result (|4.23p is a partial extension of the corresponding 
result for the zero eigenspaces (|A.35[) . Moreover, Eq. (|4.23l) implies that the 
mode spaces are covariant (see App. IA.3|) and the Lyapunov modes have well- 
defined Lyapunov exponents in the positive time limit. Hence our claimed 
tangent dynamics, via its symplecticity property (|4.20j) . implies that the 
modes are phase space perturbations which also have their (in)stability in 
the positive time limit characterized by the same Lyapunov exponent as in 
the negative limit. They therefore attain additional physical significance. 



5 Numerical Results 

5.1 Modified Numerical Scheme 

We now check the validity of the claimed tangent dynamics (|4.15p via the 
following modification of the numerical scheme. Firstly, we choose 2N lin- 
early independent tangent vectors, and allow them to converge to an or- 
thonormal set of Lyapunov vectors under the usual numerical scheme. These 
approximate Lyapunov vectors, which we denote as 6£j, are a set of positive 
Lyapunov vectors at some initial time to , and we set the negative Lyapunov 
vectors to be 

S&N-j+i = JanS^j . (5.1) 

We call d£j and S^4N-j+i for some j a symplectic pair of Lyapunov vectors. 
The ordering of the Lyapunov vectors here corresponds to the order in which 
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the numerical scheme (see Sec. 12. 3|) produces the Lyapunov vectors: most 
unstable Lyapunov vector to most stable Lyapunov vector in the positive 
time limit. This ordering is a convenient choice for the purpose of presenting 
numerical data, but differs from the ordering of the positive Lyapunov vec- 
tors 6£j, j — 1, . • • , 2N, implied in Eq. (|2.17ll . However, the ordering of the 
negative Lyapunov modes here is the same as that adopted in Eq. (|4.1|) . 

Let Q = (. .. \8£j | . . .) for j = 1, . . . , 4JV, and put the initial time to = 0. 
The middle 2K columns of Q are approximately the initial Lyapunov modes. 
We denote these modes as Srjj, with j = 2N — K + 1, . . . , 2N + K to match 
the indexing of the 6£j. 

We wish to consider the evolution of the Srjj under the normalized tangent 
dynamics only, so that the Gram-Schmidt orthonormalization, which occurs 
at regular intervals, is replaced by just normalization. Unfortunately, the 
chaotic nature of the system means that the numerics is inherently unstable 
without the Gram-Schmidt procedure, because numerical noise introduces 
components of the strongly dynamically localized Lyapunov vectors, which 
exponentially dominate the modes under further evolution. To reduce this 
numerical instability, we therefore continue to Gram-Schmidt orthonormalize 
the modes with respect to the first 2N — K Lyapunov vectors. This modified 
numerical scheme for the modes may be written explicitly as 

2N-K 

8r h {t) = c k (i m - Y, s $( s &f) D t<t> ts m ( 5 - 2 ) 

z=l 

for j — 2N — K + 1, . . . , 2N + K, where is the appropriate normalization 
constant. 

We compute the Lyapunov modes Sr)j(t) for a N — 50 particle quasi- 
one-dimensional system with periodic boundary conditions and density 0.80. 
In all the following computations we use an initial set of positive Lyapunov 
vectors 5£j, j — 1, . . . , 100, generated after 1200 collisions per particle, with 
the usual Gram-Schmidt procedure applied after each collision. 

In Fig. |4] we show the Lyapunov spectrum for this system. The zero ex- 
ponents and first three positive steps are readily apparent: the two-point 
steps consist of a two-fold degenerate Lyapunov exponent associated with 
transverse modes, whilst the single four-point step consists of a four-fold de- 
generate Lyapunov exponent associated with longitudinal modes. We label 
these steps 0, Tl, LP1, and T2 according to increasing Lyapunov exponent, 
as shown in Fig. |4j Using the Lyapunov exponent value of these steps to cal- 
culate the constants in the dispersion relation (|2.10|) , one identifies a further 
two-point step, T3, so we therefore set the number of positive modes to be 
K = 13. 

5.2 Normed Differences 

^From Eq. (|4.15|) . for index j = 88, . . . , 100 one finds that the positive Lya- 
punov modes have normalized tangent dynamics Srjj(t) = Haioo— j+i(~t)j 
while the negative modes have 

6moo-j+i(t) = S(It + J 2 oo)aioo-,+i W/[(/3i) 2 + 1]" 1/2 , (5.3) 
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Fig. 4 Partial Lyapunov spectrum for a 50 particle quasi-one-dimensional system 
with periodic boundary conditions. Lyapunov exponents Aj are indexed by j ac- 
cording to their order of calculation by the numerical scheme. Only the positive 
exponents are shown, and for the sake of the clarity of the Lyapunov steps, only the 
second fifty Lyapunov exponents are plotted. Lyapunov exponents lying on identi- 
fied Lyapunov steps are indicated by different symbols with the corresponding step 
label shown in the legend. 

where aj (t) is the jth column of A(t) . Comparison of this normalized tan- 
gent dynamics with the numerical dynamics (|4.4p reveals that the positive 
modes, 8r]j{t), should coincide with Srjj, where as usual Sij* denotes numeri- 
cal dynamics. That is, for j = 88, . . . , 100, we should have that the normed 
difference 

||<%(t)-<fyj|| 2 = 0. (5.4) 

Note that if this condition is satisfied, then the normalized tangent dynamics 
must preserve orthonormality on the positive modes 0, as predicted in Eq. 
(|4.15|) . so we obtain a measure of the orthogonality of A(t). According to Eq. 
(|4.4[) , (J2oo8rjj) t — J 2 ooS*aioo-j+i (*) and hence from Eq. (|5.3p it follows that 
we additionally expect 

||^(i) - Srf.f =2(l- m? \ l]l/2 ) (5-5) 

for the negative modes. 

FigureEK shows the norm difference \\Sr)j(t)—5rij || 2 for the positive modes. 
The zero modes' normed differences remain precisely zero, and Eq. (|5.4[) 
seems to be well-satisfied for the Tl positive modes. However the LP1 positive 
modes exhibit growth in their respective normed differences. 



2 The converse of this does not necessarily hold. 
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Fig. 5 Normed differences \\5i~ij (n) — Sr/™\\ 2 for positive (a) and negative (b) modes 
over 500 collisions, with j — 92, . . . , 97 and j = 104, . . . , 109 respectively, n is the 
collision number. The corresponding Lyapunov step label for each mode is shown 
in the legend, (a): The positive zero modes are not shown, but all satisfy \\Sr]j(n) — 

5rjj\\ = within machine precision, (b): The function f(n) = 2[1 — 1/-^/ (f3n) 2 + 1] is 
best fitted to the negative zero modes at (3 — 1/192.7 with reduced \ 2 = 3.7x 10~ 4 . 
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It is instructive to compare this growth rate with that of the non-modes. 
We therefore compute ||££j(t) — <5£j|| 2 for j — 62, ... ,67 in Fig. [6J where 
5£j(t) denotes evolution of the non-mode 8£j under the modified dynamics 
()5.2p with K = 40. The non-modes exhibit much faster, and more erratic 
growth in their normed differences than the Lyapunov modes. This suggests 
that we may consider the growth of the modes' normed differences to be 
small. Further, rotation toward the more unstable modes, due to numerical 
instability, may be responsible for the small growth in the LP1 modes' normed 
differences (see Sec. 15.31) . so a reasonable conclusion seems to be that Eq. (|5.4p 
may hold up to numerical instability. 




Fig. 6 Normed differences ||<5f,-(n) - 5Q\\ 2 for non-modes j = 62, . . . , 67 over 500 
collisions, n is the collision number. Compare this figure with Fig. [5^, noting the 
different vertical scales. 



The normed difference for the negative zero modes is in excellent agree- 
ment with Eq. (|5.5[) . as expected, and the best fit to the numerical data for 
the zero modes is obtained at /3 w 1/192.7. The negative non-zero modes' 
normed differences experience a faster growth rate than that of the zero 
modes. This is likely also due to numerical instability, but the negative mode 
normed differences nonetheless resemble the form predicted by Eq. (|5.5|) . 



5.3 Inner Products 

The form of the normalized tangent dynamics [see Eqs. (|4. 15|) and Q5.3[l ] 
means that the inner product of the symplectic pair Srjj{t) and 8r]200—j+i{t) = 
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Table 1 Mode inner products 8rji(n) ■ Sr/j(n) after n = 500 collisions, where i is 
the row index and j the column index. The column labeled (top) 88 or (bottom) 
104 contains only the single entry 1.00 and is therefore not shown. Top: positive 
mode inner products for 88 < i < j < 97. Bottom: negative mode inner products 
for 104 < i < j < 113. 



(J20oSr]j)(t) should grow as 

Svj (t) ■ S^OO-j+l (t) = [(/ ^)2 + xp/2 • ( 5 - 6 ) 

In Fig. |7|) we present the inner product growth rates for the transverse 
and longitudinal modes. The zero modes are in excellent agreement with the 
predicted functional form, and fitting obtains a time scale f3 — 1/196.3. The 
growth rate for the other modes is in good agreement with the predicted 
behavior, although mode 96 and 91 appear to exhibit a decay in their inner 
product after 800 collisions. 

The orthogonality of A(t) in Eq. (|4. 1 5[) also implies that positive modes 
should stay orthonormal to positive modes and negative modes to negative 
modes. We test this orthogonality explicitly by computing Sr/i(t) ■ 5r]j(t), for 
88 < i < j < 113. Results are shown in Table [TJ 

The off-diagonal inner products in Table [1] are mostly non-zero. In or- 
der to gain a sense of scale, it is instructive to compare the size of these 
inner products with those of the non-modes, by computing • 5£j(t) for 

61 < i < j < 69. Results are shown in Table El The inner products of the 
non-modes appear larger than that of the modes: the average absolute value 
of the off-diagonal inner product for the positive modes is 0.14 and for the 
negative modes is 0.09; however for the non-modes the average is 0.31 and 
standard deviation 0.18. Note that e.g. modes 95 and 92 both seem to be 
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Fig. 7 Inner product 5r)j(n) ■ 5ij200-j+i(n) for transverse (a) and longitudinal (b) 
modes as function of collision number n. The corresponding Lyapunov step label for 
each symplectic pair of modes is shown in the legend. For comparison the predicted 
growth rate /3n/[(/3n) 2 + l] 1 / 2 is also shown, fitted to data for the zero mode inner 
products, e.g. Srjgs(n) ■ 5r?io3(n) or 5*799(71) • <5r7io2(n), which all coincide. The best 
fit for both plots is obtained at f3 = 1/196.3 with reduced \ 2 = 4-8 x 1CT 5 . 
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rotating toward the more unstable mode 88, whilst mode 94 seems to rotate 
toward mode 90. Rotation toward more unstable modes is characteristic of 
the numerical instability in the system. Since the off-diagonal inner products 
of the modes are small in the scale of the non-mode inner products, it seems 
orthogonality of the positive (negative) modes with respect to positive (neg- 
ative) modes may be preserved by the tangent dynamics up to numerical 
instability. 
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Table 2 Non-mode inner products 5£i(n) ■ S£j(n) after n = 500 collisions, where i 
is the row index, j the column index and 61 < i < j < 69. The column labeled 61 
contains only the single entry 1.00 and is not therefore shown. 



5.4 Symplecticity Property 

It is also of interest to test the symplecticity property (|4. 20[) by computing 
the sum log(||5?7j(i)||) + log(||<5r/ 2 oo-j+i (t)\\ )• Figure [5] shows this sum for 
the modes compared to the predicted curve 1/2 \og[((3t) 2 + 1]. Once more 
the zero mode data are in excellent agreement with the predicted functional 
form (|4.20p . as expected, while the non-zero modes appear to diverge from 
it. Moreover, for the longitudinal modes 92 to 95 the sum appears to grow 
linearly, rather than logarithmically, which was a key assumption in Sec. 
(|4.4p . However, just as in Fig. [5j such divergence may also be merely due to 
numerical instability. However, one cannot completely exclude the possibility 
that the actual mode tangent dynamics differs from that presented in Eq. 
(|CT5]) . 



6 Conclusion 

We have shown that the numerical dynamics of the Lyapunov modes in the 
quasi-one-dimensional system may be represented by time-modal linear com- 
binations of an initial set of Lyapunov modes. Further, we suggested a form 
for the Lyapunov mode tangent dynamics, based on the form of the zero 
mode tangent dynamics. The inner products and normed differences of the 
Lyapunov modes predicted by this tangent dynamics resemble the numerical 
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Fig. 8 Plots of log ( 1 1 5?7j (n) 1 1 ) + log(||<577200-j+i(?i)||) as function of collision num- 
ber n for j = 92, . . . , 98. The corresponding Lyapunov step label for each mode 
symplectic pair is shown in the legend. The zero mode symplectic pairs coincide, so 
that only one (j — 98) is shown. Fitting the predicted functional form (solid line) to 
the zero mode data reveals excellent agreement: the best fit occurs at (3 = 1/190.4 
with reduced \ 2 = 3.1 x 1CT 4 . 



data, but the instability of the numerical method frustrates a more mean- 
ingful analysis. Our suggested tangent dynamics implies that the Lyapunov 
modes have well-defined Lyapunov exponents in the positive time limit. 

Acknowledgements The authors are grateful to C Angstmann and T Chung for 
discussions. 



A Theoretical Summary 

A.l Tangent Space Dynamics 

Suppose we have an r-dimensional system of N particles, with phase flow 0' : M — > 
M ~ R 2rJv . Classification of the stability of a tangent vector 5£ £ T^M is achieved 
by following theorem originally due to Oseledec [?,?]. 

Theorem 1 (Multiplicative Ergodic Theorem) Let cj> : M —> M be a (not neces- 
sarily invertible) differentiate phase flow, which preserves the Liouville measure. 
Then for almost all £ € M (with respect to this measure) the stability matrix 

At,+ = hm [{D^f D^] 1/24 (A.l) 

exists and there exists a positive integer I and distinct real numbers A^' such that 
A^ !+ has eigenvalues exp(A' J ') for j = 1, . . . ,1, where > . . . > Let J7^| be 
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the eigenspace of A^ t + corresponding to eigenvalue exp(A"') 3 and let 

V ( U l = Ugl (A.2) 
y_(i+i)|Qj Then furthermore, 

V£ - = e T e M : hm \ log = A« } . (A.3) 

As previously, the numbers A"' are called Lyapunov exponents, and the ordered set 
{A (J) } is called the Ly apunov spectrum. We have included the denominator in the 
logarithm of Eq. (|A.3[) to indicate explicitly that the Lyapunov exponent describes 
a ratio of tangent vector lengths. This denominator contributes a zero term to the 
exponent in the limit t — > oo, so we neglect it henceforth. 

The stability matrix A^ t+ is clearly symmetric, so the eigenspaces U^l form an 
orthogonal decomposition of the tangent space. That is, T^M — UjQ © • • • © 

and (S£,5r]} — for 8£ £ ^i+i ^ ^e, k + anc ^ 3 ^ ^, wnere we use the usual 
inner product here - the dot product - which defines the norm || ■ ||. It follows from 
this decomposition and the definition of the spaces Vf?\ that T^M — D ... D 
v m D yil+l) = {0} . 

According to the theorem, the sets VPJ. — + ^ partition T^M according to 
the Lyapunov exponents: each non-zero tangent vector has an unique associated 
Lyapunov exponent. This means that we may write 

k—j 

= \8£,eTcM: lim i log IIDf^ll < A (j) lu{0} , (A.4) 

[ t— >oo t J 

and call the (I — j + l)-th most stable subspace of T(M in the positive time 
limit, or the positive (I — j + l)-th stable subspace for short. From the definition 
of the spaces in the theorem it also follows immediately that 

and hence the non-zero elements of the eigenspaces have associated Lyapunov ex- 
pone nt A "', as would be expected. Further, orthogonality of the eigenspaces and 
Eq. ()A.2[) implies 

uft. = v& n (V^Y , (A.6) 

so that the eigenspace E/j^. is that subspace orthogonal to the positive (I — j)-th 
stable subspace but contained within the positive (I — j + l)-th stable subspace. 

A.2 Hamiltonian Properties 

The phase flow of the particle system is Hamiltonian, and consequently the com- 
position of phase flows exhibits the (homomorphism) property 



0* o (f) s = 0* 



(A.7) 
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The phase flow is clearly invertible, with inverse (</>*) 1 = cj> ' an d cj>° is the identity 
map. (As such, the set {<^*}tgK is a n Ab elian group that is isomorphic to the real 
numbers.) We have then from Eq. (IA.7|l and the chain rule that the tangent flow 
is a co-cycle, that is 

D^s^Drf* = D ( d> t+S . (A.8) 
The tangent flow is therefore invertible, with inverse (D^rjf)^ 1 = D^t(f-\<f>~ t and 
D^(f>° — I, the identity matrix. 

A Hamiltonian phase flow preserves the Liouville measure on M [?] , so that the 
phase flow under consideration satifies the hypotheses of the Multiplicative Ergodic 
Theorem. Invertibility of the phase flow means that we may apply the theorem also 
to the inverse <j> , which is equivalent to contemplating a switch in the direction 
of time. It follows from the theorem that the stability matrix 

Ae-= lim [(D 5 t ) T De0 t l 1/2|t| (A.9) 

t— y— oo L J 

also exists, and A^,- must have the inverse eigenvalues to , so that it has 
eigenvalues exp(— A' J ^), for j = 1, . . . , I. As we expect, the Lyapunov exponents of 
_t are then the opposite sign to those of $ but have the same multiplicity, denoted 
mult(j). We let U^l be the eigenspace corresponding to eigenvalue exp(— A"'), and 
then T^M = ?7j! © • • ■ © Uc _ is another orthogonal decomposition of the tangent 
space. In general, note that / Jji _ , but we do have that dim J/y| = mult(j) = 
dim U^_. 

The j-th most stable subspace in the negative time limit VJr _ , called the nega- 
tive j-th stable space, is now well-defined by the theorem. Since — A' 1 ' < . . . < — A''' 
are the Lyapunov exponents, we must hav«Q 

Vgl = ujV®'-QUgl , (A. 10) 

and then the tangent space is partitioned into the sets 

V U)_ _ V U-V = |^ G T £ M : t lim -i- log [[D^S^W = -A (j) | , (A.ll) 
where T e M = D ... D V^l D V^l = {0}. From Eq. jX30) we have 

jjW = V U)_ n (y^Y , (A.12) 

similarly to Eq. ()A.6|I . The non-zero elements of the eigenspaces Ug _ has associated 

exponent —A"' analogously to Eq. (|A.5[) , and note also that, just as in Eq. (|A.4[I . 
here 

V e u l = e T e M : lim A lo s llAr^ll < u {o} . (A.13) 

I 1*1 ' J 

We call the elements of the eigenspaces 1/1 + (U^}_) the Lyapunov vectors, each 
of which has an associated Lyapunov exponent in the positive (negative) time limit. 
We may clearly choose an orthonormal basis of T^M such that each element of the 

basis is an element of an eigenspace £/! + {U^l.), just as claimed in Section [2.21 
Further, if any of the Lyapunov exponents have multiplicity greater than unity, then 
we have many choices of such a basis and only the eigenspace itself is well-defined. 

3 Note that since the order of the Lyapunov exponents has been reversed, we 
have reversed the indexing of these spaces so that the + l)-th stable space is 

y(l-o+n rather than y(ri_ 
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A. 3 Covariant Subspaces 

Suppose the tangent vector <5£ £ Vj + — V}^ 1 ^, so that it has corresponding Lya- 
punov exponent A' J '. If r\ = S (£) then 5rj = D^^S^ G T^M should also have 
Lyapunov exponent since the Lyapunov exponents are independent of time. 
This property is guaranteed by the co-cycle property of the tangent flow and its 
consequent invertibility. It can be shown that 

V$ &>+ = DtfVgl , (A.14) 

for s finite. Further, since the V^l are closed subspaces, the equality also holds in 
either limit s — > ±00. A subspace satisfying Eq. (|A,14[) for all real s and s — ► ±00 
is called a covariant subspace, and an identical result holds for the . Note that 

the set V^J_ — V^ 1-1 ' is not covariant since the equality does not necessarily hold 
in the limit s — » —00, as this set is not closed. 
For a non-zero 

<5£ G V^+, we have from Eq. jAjgl that D^- s 5^ £ V^ s( + . 

Substituting this into the inequality in Eq. (|A.4[l and applying the successive trans- 
formations t — * t — s then s — *• — (s + t), which together are the map t — > — s, one 
finds that 



V® C \S£ e T t M : t Ita log WD^S^W > -A (j) |. U {0} , (A.15) 



and similarly 



VP c\d£& TtM : lim -!- log W^S^ > X U) } U {0} . (A.16) 
5 ' [ |i| J 

Equation (|A,5[) implies the eigenspace Ujr*^_ (U^}_) has associated Lyapunov 

expo nent (— A' J ') in the positive (negative) time limit. However, Eqs. (|A.15|) 
and (|A.16|) unfortunately imply that in the negative (positive) time limit we can 

only say that a non-zero element of (U^}_) has exponent A > — A*- 3 ' (A > A' J '). 
Nevertheless, comparison of Eqs. ()A.4[) . (|A.13|) . (|A.15() and (|A.16|) suggests that 
the intersection of the spaces and Ve_ produces a subspace with well-defined 
Lyapunov exponent in both the positive or negative time limit. That is, let 



Then 



W U) _ V U) n yU) (A 17) 



W^ = {8teT e M: t ]^±] a g\\D e <l, t 8Z\\=±xM) j{oj . (A.Js, 



and it can be proven (see e.g. Ref. [?]) that 

vgl = wf 5 • • • W^ j) , (A. 19) 

so that T 5 M = W^ 1 ' • • ■ W* W . This is called the Oseledec decomposition of the 

(AS 

tangent space, and W£ is called an Oseledec space. Note that the Oseledec spaces 
are not orthogonal like the but are clearly covariant spaces unlike the U^_. 

From Eqs. (|AT9)l and (|AT2)l it follows that dim Ui j) _ = dim W { J ] . 
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A. 4 Numerical Scheme 

The Lyapunov exponents are calculated via the numerical scheme of Benettin [?] 
and Shimada [?]. A byproduct of this scheme is a orthonormal basis of Lyapunov 

vectors that are elements of the eigenspaces Uir_ . The numerical scheme depends 
upon the following well-known result, which we state as a proposition. We call the 
map Z?| 0' 1| Z?£ || the normalized tangent flow. 

Proposition A randomly chosen non-zero, tangent vector converges for t — > +oo 
under the normalized tangent flow to the most stable space in the negative time 

limit, That is for non-zero <5£ G T^M and t —t oo, 



¥wm\ ' proj r^)-pS%J' (A - 20) 

in the norm || • ||, where proj(V, x) denotes the projection of the vector x onto the 
subspace V. 

The proposition can be proved by use of the Oseledec decomposition of the 
tangent space, the invariance property of the Oseledec spaces and the fact that 
these spaces possess well-defined, ordered Lyapunov exponents (see e.g. [?]). The 
idea here is that a randomly chosen tangent vector converges to its projection onto 

the space = under the normalized tangent flow, but it does not converge 

at all in the norm || • | under the usual tangent flow. 
Let 

M(j) = j k e Z : mult(i) < k < ^ mult(i) I . (A.21) 

1=1 1=1 ' 

The proposition may be generalized to show that for k randomly chosen, nonzero, 
linearly independent tangent vectors 5£i, 

c ^ span {p^&}L Q ■ (A ' 22) 

where j is such that k £ M(j). 

We now fix t and define {Sr/"} recursively as the Gram-Schmidt orthonormalized 
set generated by {[-D^n-iit^)^*]^"^ 1 }, where {dr/f} is the orthonormalized set of 
the tangent vect or set Note that the tangent flow is not orthogonal. It follows 

then from Eqs. (|A.10|) and (|A.22|) that 

lim ^ e lim E$> _ (A.23) 

for k £ N(j). In other words, the numerical scheme operates as follows: if we 
randomly choose a set of 2rN linearly independent tangent vectors, then apply 
the tangent flow to these vectors whilst regularly Gram-Schmidt orthonormalizing 
them, we produce in the positive time limit an orthonormal basis of Lyapunov 
vectors. These Lyapunov vectors have well-defined Lyapunov exponents only in the 
negative time limit. 



A. 5 Further Hamiltoman Properties 

A Hamiltonian phase flow also preserves the symplectic structure of the phase space 
M, so that the tangent flow is a symplectic matrix. Hence 

(D^f J2rND ( c/> i = J2rN , (A.24) 
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where J2 r jv is defined in Eq. (I4.2[) . and we put J = J2rN henceforth. It is well- 
known that, due to this symplecticity, there is always a zero Lyapunov exponent 
and all other exponents are paired signwise with equal multiplicity such that 

X {i-i+V) = _ A 0) i (A. 25) 

Hence the number of distinct exponents, I, is always odd and ^(I'+ 1 1/ 2 ) = 0. Hence- 
forth we let m = (I + l)/2, and we call the associated eigenspaces U^} or U^} 

the zero eigenspaces. 

The stability matrices Af,± (|A.1IA.9|I are thus also symplectic. It can then be 
shown from Eq. ()A.25[) and the definition of the as eigenspaces of that 

U\% = JU^ +1) (A.26) 
and hence for the zero eigenspaces 

u t2 = JU< C± ■ ( A - 27 ) 

Note that for all S£ G T e M, the dot product <5£ ■ J8£ = 0. We deduce dim = 
mult(m) must be even, since if <5£ is an element of an orthogonal basis for U^±, 

then by Eq. (fA~27|l so is JS£. 

An important consequence of Eqs. (|A.25[) and (|A.26|l is that the Lyapunov 
spectrum and an orthonormal basis of Lyapunov vectors is fully described by con- 
sidering only positive Lyapunov exponents and the Lyapunov vectors corresponding 
to these. That is, for a orthonormal set of Lyapunov vectors {5^}, k = 1, . . . , 2rN, 
we may choose S&rN-j+i = J8£j for j = 1, ... , rN. 



A. 6 Symmetries and Zero Modes 

If we presume that the dynamics of the system under consideration has full (trans- 
lational) spatial and time symmetry, then there are consequently r + 1 tangent 
vectors <5£i, i = 1, . . . , r + 1 such that 

DztfSZi = 5& . (A.28) 

The first r tangent vectors of these correspond to a uniform tra nslat ion of the 
N particles in the i-th spatial direction. With reference to Eq. (|2.3p these are 
8£i = l/TV 1 / 2 (0, d, ■ ■ . , ei) , where G R r is the i-th unit basis vector repeated 
N times and is the rN dimensional zero vector. Symmetry in time means that 
the system also possesses time translational invariance, and hence the velocity unit 
tangent vector at state £ = (p, q), 8£ r+ i = l/||p|| (0,p), satisfies D^tftSt^r+i = 5£ r +i- 
By Noether's theorem [?], the spatial and time symmetries of the Hamiltonian 
dynamics give rise to conserved quantities, namely momentum and energy respec- 
tively. Hence we have a further r tangent vectors corresponding to a momentum 
shift in the i-th direction <5£i+ r +i = 1/N 1 ^ 2 (e,, . ..,6i, 0), for i = 1, . . . ,r. These 
tangent vectors have dynamics 

£»e0%+r+i = Sii+r+l + Pt5& , (A.29) 

where j3 > is some constant, such that they are generalized eigenvectors of D^c/y 
with algebraic multiplicity 2. Further, a shift in energy effectively raises or lowers 
the temperature of the system, and hence we have a tangent vector <5^2r+2 = 
(p, 0) which has dynamics 



(A.30) 
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corresponding to linear separation in time. Note that we have written all the 5£i 
with unit norm. 

Now, 5£i+ r +i = —JS^i for all i = 1, . . . , r + 1. Applyi ng this relatio n toge ther 
with the symplecticity of the tangent flow to Eqs. (IA.28j) . (|A.29j) and (|A.30[) one 
finds that for i = 1, . . . , r + 1 

(D^fS^ = 5& + J/3t5& , (A.31) 

ft follows by the def inition of t he sta bility matrices in the Multiplicative Ergodic 
Theorem with Eqs. (|A.28|I and (|A.31f) that 

I", nT *T 1/21*1 

4e,±tf&= lim fl^"i) e fl 5^ 

i— >-±oo L J 

= lim (I + f3tj) 1/m S& 

t— »±oo 

= % , (A.32) 

for i = 1, . . . , r + 1, and a similar result holds for i = r + 2, . . . , 2r + 2. By definition 
of [/i™ and l/i _ as a zero eigenspace of and vl^- respectively, we have 

span{^} 2 ^ 2 C n U$ . (A.33) 
We therefore call S£i the zero modes. 



A. 7 Equivalence of Zero Spaces 

Empirical results confirm that the linear combinations of the S(,i i n Ap p. IA.6I are 
indeed the zero modes generated by the numerical scheme of Add. |A~4| Moreover, 
one finds the multiplicity of A' m \ mult(m) = 2s + 2, for s the number of unbroken 
spatial symmetries in the system of interest. It follows that dim L^™' = dim Uc_ = 
dim(sp an{(5gi }). Since the zero modes are linearly independent, then with reference 
to Eq. ()Al33[) it must be that 

U^J = span{^} 2r+2 = tfi+ J , (A.34) 

In other words, the zero eigenspaces coincide and their non-zero elements conse- 
quently have a well-defined Lyapunov exponent A = in either the positive or 
negative time limit. 

^From the definition of the subspace VPJ. in the Multiplicative Ergodic Theo- 
rem, Eq. CPU implies that llg} C V^ '. Hence f/£_° C V^™ } n v£_° = W\ m) . 
However, dim = dim so it must be that 

w ( m ) = ^ = p<») ( ( A .3 5 ) 

or in other words the zero eigenspaces and the zero Oseledec subspace coincide. 
This result is a direct consequence of the structure of the zero modes and the 
empirically verified multiplicity of the zero exponent. 
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B Symplecticity of Numerical Dynamics 

Consider the numerical dynamics fi 1 = Nq{0) denned in Sec. 13.11 and let Q — 
(. . . \5idj\ . . .) be an orthogonal 2n x 2n matrix of Lyapunov vectors such that 
5uj2n-j+\ = JSuj. Note that this choice of Sujj coincides with a possible choice 
of an ordered set of Lyapunov vectors generated by the numerical scheme (see 
App. IA.5|I . We now show that Nq is a symplectic operator, in the sense that 

JNhiSwj) =Nh(J6u s ) , (B.l) 

for j = 1, . . . , In and J = J2 n as defined in Eq. (|4.2|l . 

Firstly, we may consider the numerical dynamics to consist of evolution under 
the tangent flow followed by the Gram-Schmidt orthonormalization proce- 

dure(see Add IA.4"]) . Letting Q(t) = D^Q, then the numerical dynamics may be 
written as £2 f — J?(f)i? _1 (f) where R(t) is the upper triangular matrix of the so- 
called QR factorization of S7(t). Hence we have 

Nh{n) = D i <t> t C2R- 1 {t) , (B.2) 

as a decomposition of the numerical dynamics. A"^ is explicitly dependent on the 
choice of Q and is implicitly dependent on the ordering of the 8u)j due to the 
iterative nature of the Gram-Schmidt procedure, so it is not a linear map, although 
the matrix Nq{Q) is clearly orthogonal. 

Now, (Sujf = Na(5ujj) = D i 4> t flR~ 1 (t)fl T 5uj J . Defining the matrix 

N(t) = D e <f> t nR~ 1 {t)n T , (B.3) 

we also have that N n (JSuij) = Nn(8 m2n-j+i) ~ N(t)J5ujj. Since the 5ujj form a 
basis of the tangent space, Eq. (|B.1|) is then equivalent to JN{t) — N(t)J. Since 
N(t) is or thogo nal, it therefore suffices to show that N(t) is symplectic in order to 
show Eq. (IB. If) . 

Let TZ(t) = nR{t)Q T . Then the orthogonality of N(t), i.e. N^it) = N T (t), 
together with Eq. (|Fl3)) implies that Tl(t) T 1Z{t) = (D^) 7 ' . It follows that 
1Z(t) T 1Z(t) is symplectic due to the symplecticity of the tangent flow. That is, 

H(t) T H(t)J1l(t) T 1l(t) = J . (B.4) 
By hypothesis, we may write 

n = (zi|jzi/) (B.5) 

where A is an In x n matrix with orthonormal columns that satisfies A T J A — 0, 
and I is the n x n antidiagonal matrix with Is on the antidiagonal and zeros 
elsewhere. It then follows from this choice that 

Q T JH = J , where J = (j • (B.6) 

Hence Eq. (IB.4[) implies 

R T (t)R{t)JR T {t)R(t) = J . (B.7) 

By E q. (IB.3|) . in order for N(t) to be symplectic we require TZ T JTZ = J. From 
Eq. (|B.6|) this requirement is equivalent to 



R(t) T JR{t) = J. 



(B.8) 
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It is therefore s uffici ent to show Eq. (|B.8[l in order to show that N(t) is symplectic 
and hence Eq. (|B.ip holds. 

One can show via some algebr a tha t Eq. (|B.T|) together with the upper trian- 
gular nature of R(t) produces Eq. (|B,8[l . In detail, we write 

R(t) = (B.9) 

where A, C are themselves n x n upper-triangular matrices with strictly positive 
values on the diagonal. Eq. (|B.7|I provides restrictions on A, B and C so that one 
finds 

R{t)JR T {t) = ( C/ ^ T ~ AI ° T ) ■ (B.10) 

AIC T is an upper antitriangular matrix (an upper triangular matrix multiplied 
by I) with strictly positive values on the antidiagonal,and one finds further that 
AIC T is also orthogonal. Th e onl y possibility is AIC T = I, so since J 2 = —7, Eq. 
(|B.8|) holds and so does Eq. (|B.1|> . 



